R Notebook

For generating most, but not all figures in the manuscript.

library(tidyverse)
library(cowplot)
library(lubridate)
library(mgcv)
source("UVP_2017_library.R")

Parsed with column specification:
cols(
  Cruise = col_character(),
  Station = col_character(),
  `mon/dd/yyyy` = col_character(),
  `hh:mm` = col_time(format = ""),
  `Longitude [degrees east]` = col_double(),
  `Latitude [degrees north]` = col_double(),
  `Bottom Depth [m]` = col_double(),
  `Pressure [db]` = col_double(),
  `Temperature [degrees C]` = col_double(),
  `Temperature 2 [degrees C]` = col_double(),
  `Salinity [psu]` = col_double(),
  `Salinity 2 [psu]` = col_double(),
  `Fluorescense [mg/m^3]` = col_double(),
  `Beam Transmission [%]` = col_double(),
  PAR = col_double(),
  `Oxygen [umol/kg]` = col_double(),
  `Oxygen [% saturation]` = col_double()
)

theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')

Particles Only

Read In Data

These read back in in UTC, since write_csv saves everthing into that time zone.

options(readr.default_locale=readr::locale(tz="Mexico/General"))
bes<- read_csv("dataOut/binned_EachSize.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.

bds <- read_csv("dataOut/binned_DepthSummary.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.

ues <- read_csv("dataOut/unbinned_EachSize.csv")

Parsed with column specification:
cols(
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = ""),
  depth = col_double(),
  psd_gam = col_double(),
  vol = col_double(),
  sizeclass = col_character(),
  lb = col_double(),
  ub = col_double(),
  binsize = col_double(),
  TotalParticles = col_double(),
  nparticles = col_double(),
  n_nparticles = col_double(),
  biovolume = col_double(),
  speed = col_double(),
  flux = col_double(),
  flux_fit = col_double(),
  GamPredictTP = col_double()
)

uds <- read_csv("dataOut/unbinned_DepthSummary.csv")

Parsed with column specification:
cols(
  .default = col_double(),
  project = col_character(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.

Specify the base of the photic zone (which is inside of the OMZ) and the bae of the OMZ

PhoticBase <- 160
OMZBase <- 900
DVMBase <- 600
OMZTop <- 90

Figure 4

Total Particle numbers and particle size distribution slope

library(scales)
#https://stackoverflow.com/questions/30179442/plotting-minor-breaks-on-a-log-scale-with-ggplot
log_breaks = function(maj, radix=10) {
  function(x) {
    minx         = floor(min(logb(x,radix), na.rm=T)) - 1
    maxx         = ceiling(max(logb(x,radix), na.rm=T)) + 1
    n_major      = maxx - minx + 1
    major_breaks = seq(minx, maxx, by=1)
    if (maj) {
      breaks = major_breaks
    } else {
      steps = logb(1:(radix-1),radix)
      breaks = rep(steps, times=n_major) +
               rep(major_breaks, each=radix-1)
    }
    radix^breaks
  }
}
scale_x_log_eng = function(..., radix=10) {
  scale_x_continuous(...,
                     trans=log_trans(radix),
                     breaks=log_breaks(TRUE, radix),
                     minor_breaks=log_breaks(FALSE, radix))
}

ybreaks = seq(from = 1200, to = 0, by = -50)
ylabels = ybreaks
ylabels[c(FALSE, TRUE, TRUE, TRUE)] <- ""

#theme_set(theme_bw)


#theme_set(theme_cowplot)

PlotParticlesmany <- uds %>% 
  filter(project == "ETNP") %>%
  ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
 
  #geom_path(aes(x = psd_gam)) + 
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
  geom_point(alpha = .6, size = 2, stroke = 1) +
  scale_y_reverse(limits = c(1200, 0),
                  breaks = ybreaks,
                  labels = ylabels) +
  scale_shape_manual(values = c(21:25)) +
  scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black"), limits = c(0, 24)) +
  scale_x_log10(breaks = c(10, 100, 1000), minor = c(5, 50, 500)) +
  #theme(legend.position = "none") +
  #scale_x_log_eng()+
  labs(y = "Depth (m)", x = "Particles / L") + 
  geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = OMZBase, color = "darkblue") +
  geom_hline(yintercept = OMZTop, color = "darkblue") +
  #geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed") +
  theme(legend.position = "none")

PlotPSDmany <- uds %>% 
  filter(project == "ETNP") %>%
  ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
 
  #geom_path(aes(x = psd_gam)) + 
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
  geom_point(alpha = .6, size = 2, stroke = 1) +
  scale_y_reverse(limits = c(1200, 0), breaks = ybreaks, labels = ylabels) +
  scale_shape_manual(name = "Day", values = c(21:25)) +
  scale_fill_gradientn(name = "Hour",breaks = c(0, 6, 12, 18, 24), labels =  c("0", "6", "12", "18", "24"), colors = c("black", "blue", "white", "orange", "black"), limits = c(0, 24)) +
  labs(y = "Depth (m)", x = "Particle Size Distribution Slope") + 
  geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = OMZBase, color = "darkblue") +
  geom_hline(yintercept = OMZTop, color = "darkblue") +
  #geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed") +
  coord_cartesian(clip = "off") +
  annotate(geom = "text", x = -2, y = 1100, label = "More large \n particles", inherit.aes = FALSE, size = 6) +
  annotate(geom = "text", x = -4.5, y = 1100, label = "More small \n particles", inherit.aes = FALSE, size = 6)

Ignoring unknown parameters: inherit.aesIgnoring unknown parameters: inherit.aes

  # annotation_custom(
  #     grob = textGrob(label = "More big particles", hjust = 0, gp = gpar(cex = 12)),
  #     ymin = 0,      # Vertical position of the textGrob
  #     ymax = 0,
  #     xmin = 0,         # Note: The grobs are positioned outside the plot area
  #     xmax = 0)
  

plot_grid(
  PlotParticlesmany,
  PlotPSDmany,
  rel_widths = c(1, 2),
  labels = c("A", "B")
  )

Removed 266 rows containing missing values (geom_point).Removed 266 rows containing missing values (geom_point).

ggsave("figures/ParticlesPSDMany.png")

Saving 12 x 7.41 in image

ggsave("figures/ParticlesPSDMany.svg")

Saving 12 x 7.41 in image

Figure 3 Summary Statistics

Particle number vs depth and time

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

#FSG4 <- gam(tot_nparticles~ s(depth, k = 3)  + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

summary(FSG1)

Family: gaussian 
Link function: identity 

Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, 
    bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.96538    0.09655   92.86   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
            edf Ref.df     F p-value  
s(depth) 1.0000  1.000 5.663  0.0202 *
s(Day)   1.4473  1.694 2.033  0.0921 .
s(Hour)  0.6402  2.000 0.502  0.2154  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.135   Deviance explained = 17.4%
GCV = 0.68369  Scale est. = 0.64319   n = 69

summary(FSG2)

Family: gaussian 
Link function: identity 

Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  8.96538    0.09727   92.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value  
s(depth) 1.000  1.000 5.741  0.0194 *
s(Day)   1.469  1.718 2.214  0.0792 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.122   Deviance explained = 15.3%
GCV = 0.68744  Scale est. = 0.65288   n = 69

summary(FSG3)

Family: gaussian 
Link function: identity 

Formula:
tot_nparticles ~ s(depth, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   8.9654     0.1004   89.34   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
         edf Ref.df     F p-value  
s(depth)   1      1 5.736  0.0194 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.0651   Deviance explained = 7.89%
GCV = 0.71559  Scale est. = 0.69485   n = 69

#summary(FSG4)

summary(FSG1)$r.sq - summary(FSG2)$r.sq # extra R^2 explained by hour

[1] 0.01304448

summary(FSG2)$r.sq - summary(FSG3)$r.sq # extra explained by day

[1] 0.05646094

summary(FSG3)$r.sq # R^2 explained by depth

[1] 0.06511156

There is no statisticlly significant affect of time on particle number. If I take all of the time variables out and just compare to depth, there is a relationship to depth p = 0.02. But the R^2 is only 6.5%. Pretty weak

Particle size distribution vs depth and time

bdsAddTime <- bds %>%
  mutate(Hour = hour(time), Day = day(time))

FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

FSG4 <- gam(psd~ s(depth, k = 3)  + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

summary(FSG1)

Family: gaussian 
Link function: identity 

Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.90483    0.01809  -215.8   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F p-value    
s(depth) 1.848  1.977 88.789  <2e-16 ***
s(Day)   1.506  1.756  0.953  0.4853    
s(Hour)  1.361  2.000  2.260  0.0424 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.716   Deviance explained = 73.6%
GCV = 0.024631  Scale est. = 0.022591  n = 69

#summary(FSG2)
#summary(FSG3)
summary(FSG4)

Family: gaussian 
Link function: identity 

Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.90483    0.01814  -215.2   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F p-value    
s(depth) 1.845  1.976 88.630  <2e-16 ***
s(Hour)  1.524  2.000  2.559  0.0409 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.714   Deviance explained = 72.8%
GCV = 0.024244  Scale est. = 0.022709  n = 69

summary(FSG3)$r.sq # R^2 of depth

[1] 0.6929817

summary(FSG4)$r.sq - summary(FSG3)$r.sq # Improvement from adding hour of the day

[1] 0.02133449

summary(FSG1)$r.sq - summary(FSG4)$r.sq # Improvement from then adding day of the week

[1] 0.001481684

PSD varies with depth, but doesn’t statistically relate to hor orday. Comparing the R2 values from models tells us that you explain 69% of varience with depth.

Figure S6

Comparing the two stations

S6A Number vs depth

PlotNParticlesEP <- uds %>% 
  filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
 geom_point(alpha = 0.7, size = 2, stroke = 1) +
  #geom_path(aes(x = tot_nparticles)) +
  #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
  labs(x = "Particles/L", y = "Depth (m)") +
  theme(legend.position = "none") +
  scale_shape_manual(values = c(1:5)) +
  geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = 200, color = "darkgreen") +
  geom_hline(yintercept = OMZBase, color = "darkblue") 

PlotNParticlesEP

I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)

S6B Particle size distribution vs depth

PlotPSDEP <- uds %>% 
  filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
 geom_point(alpha = 0.7, size = 2, stroke = 1) +
  geom_path(aes(x = psd_gam)) +
  geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_color_manual(values = c("gray20", "brown"))  +
  scale_shape_manual(values = c(1:5)) + labs(y = "", x = "Particle Size Distribution Slope") +
  geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = 200, color = "darkgreen") +
  geom_hline(yintercept = OMZBase, color = "darkblue") 

PlotPSDEP

Figure S6 Combined

plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))

Removed 1211 rows containing missing values (geom_point).Removed 1211 rows containing missing values (geom_point).Removed 1211 row(s) containing missing values (geom_path).

ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")

Saving 10 x 4 in image

ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")

Saving 10 x 4 in image

Figure S7

Large and spall particle number, flux and size

mainParticleComponents <- bds %>%
  filter(profile %in% c("stn_043", "p16n_100")) %>%
  select(project, profile, depth,
         tot_nparticles, small_nparticles, big_nparticles,
         tot_psd = psd, small_psd, big_psd,
         tot_flux_fit, small_flux_fit, big_flux_fit) %>%
  pivot_longer(cols = -c("project", "profile", "depth")) %>%
  separate(name, c("size", "meas")) %>%
  mutate(meas = recode(meas, nparticles = "particles/L")) %>%
  mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))

Expected 2 pieces. Additional pieces discarded in 273 rows [7, 8, 9, 16, 17, 18, 25, 26, 27, 34, 35, 36, 43, 44, 45, 52, 53, 54, 61, 62, ...].

PlotFlx <- mainParticleComponents %>% 
  filter(meas != "psd") %>%
  ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size ~ meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + theme(axis.title.x = element_blank(), legend.position = "none", strip.background.y = element_blank(), strip.text.y = element_blank(), plot.margin = unit(c(7,0,7,7), "pt")) + scale_color_manual(values = c(P16 = "brown", ETNP = "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = PhoticBase, color = "darkgreen")

PlotPSD <- mainParticleComponents %>% 
  filter(meas == "psd") %>%
  ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size~meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = unit(c(7,7,26.5,0), "pt")) +
  scale_color_manual(values = c(P16 = "brown", ETNP = "gray20")) +  scale_shape_manual(values = c(1:5)) +  theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = PhoticBase, color = "darkgreen")

plot_grid(PlotFlx, PlotPSD, rel_widths = c(3, 2))

Removed 246 rows containing missing values (geom_point).Removed 123 rows containing missing values (geom_point).

ggsave("figures/BigVsSmall.svg")

Saving 7.29 x 4.5 in image

ggsave("figures/BigVsSmall.png")

Saving 7.29 x 4.5 in image

Flux small and flux tot track so closely because particle fractal dimension alpha, plus flux fractal dimension, gamma > |psd|. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top.

Figure S4

Example particle size distributions

eg_dataline <- bds %>% 
  filter(profile == "stn_043", depth == 162.5)
eg_slope =  eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)

eg_datablock <- bes %>%
  filter(profile == "stn_043", depth == 162.5)


eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)

eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)


EgNNP <- eg_datablock %>%
  ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() + 
  geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)", x = "Size (mm)")

EgNP <- eg_datablock %>%
  ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() + 
  geom_path(data = eg_df) + labs(y = "Normalized Particles" , x = "Size (mm)")

EgTP <- eg_datablock %>%
  ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() + 
  geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)", x = "Size (mm)")

plot_grid(EgNNP, EgTP, labels = c("A", "B"))

Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis

ggsave("figures/ExamplePSD163m.png")

Saving 7.29 x 4.5 in image

ggsave("figures/ExamplePSD163m.svg")

Saving 7.29 x 4.5 in image

Figure 5A

Flux attenuation with respect ot depth and time. All extrapolated from the UVP and traps combined.

scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_b <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }

scientific_10_c <- function(x) {
    xout <- gsub("1e", "10^{", format(x),fixed=TRUE)
    xout <- gsub("{-0", "{-", xout,fixed=TRUE)
    xout <- gsub("{+", "{", xout,fixed=TRUE)
    xout <- gsub("{0", "{", xout,fixed=TRUE)
    xout <- paste(xout,"}",sep="")
    return(parse(text=xout))
    
}

scale_x_log10nice <- function(name=NULL,omag=seq(-10,20),...) {
    breaks10 <- 10^omag
    scale_x_log10(breaks=breaks10,labels=scientific_10_c(breaks10),...)
}


#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time)))  + geom_point(size = 2, stroke = 1)+
  #geom_path() +
  scale_y_reverse(limits = c(1000, 0), breaks = ybreaks, labels = ylabels) +
  scale_x_log10nice()+
  #scale_x_log10() + 
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + 
  scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black"), limits = c(0, 24)) +
  
  labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
  #labs(x = "moo", y = "Depth (m)") +
  geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE) +
  theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm")) +
   geom_segment(aes(y = PhoticBase, yend = PhoticBase, x = 20, xend = 500), color = "darkgreen", stroke = 0.5)+
   geom_segment(aes(y = OMZBase, yend = OMZBase, x = 20, xend = 500), color = "darkblue", stroke = 0.5) +
   geom_segment(aes(y = OMZTop, yend = OMZTop, x = 20, xend = 500), color = "darkblue", stroke = 0.5)

Ignoring unknown parameters: strokeIgnoring unknown parameters: strokeIgnoring unknown parameters: stroke

  #geom_segment(aes(y = DVMBase, yend = DVMBase, x = 20, xend = 500), color = "darkgoldenrod", lty = "dashed", stroke = 0.5)
  
  #+ geom_hline(yintercept = OMZBase, color = "darkblue")



pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlxLegend <- get_legend(pltFlx)

Removed 14 rows containing missing values (geom_point).

pltFlx

#plotly::ggplotly(plt1)

Figure 5B

Zooming in on where the action is happening

ylabels_fine <- ybreaks
ylabels_fine[c(FALSE, TRUE)] <- ""
#xbreaks_fine <- c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20))
xbreaks_fine <- seq(from = 20, to = 180, by = 10)
xlabels_fine <- xbreaks_fine
toss <- c(rep(FALSE, 5), rep(c(TRUE, FALSE), 6))
xlabels_fine[toss] <- ""
xlabels_fine[xlabels_fine %in% c(120, 160)] <- ""

#lenx <- length(xlabels_fine)
#xlabels_fine[c(lenx, lenx-2)] <- ""
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
  #geom_path() +
  scale_y_reverse(limits = c(500, 75), breaks = ybreaks, labels = ylabels_fine) +
  #scale_x_log10() +
  scale_x_log10(breaks = xbreaks_fine,
                labels = xlabels_fine,
                limits = c(20, 180)) +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) + 
  scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux", y = "Depth") + theme(legend.position = "none")+
geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = OMZTop, color = "darkblue") +
  theme(axis.text.x = element_text(vjust = .5))

pltFlxZoom

#plotly::ggplotly(plt1)

Figure 5C

Rate of change of flux, taken to the fifth root so one can see patterns.

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time)))  + geom_point(size = 2, stroke = 1)+
  #geom_path() +
  scale_y_reverse(limits = c(1000, 0), breaks = ybreaks, labels = ylabels) +
  scale_x_continuous(limits = c(-2.1, .6), breaks = seq(from = -2, to = .75, by = 0.5)) +
  #scale_x_log10() +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + 
  scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  geom_vline(xintercept = 0) +
  labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")+
  geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = OMZBase, color = "darkblue") +
  geom_hline(yintercept = OMZTop, color = "darkblue")
  #geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed")
  #labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")

pltDelta3

#plotly::ggplotly(plt1pos)

Combining the plots

Within panel drawing

pgTop <- ggdraw(pltFlxNoLegend 
       ) +
  draw_plot(pltFlxZoom, .4, .25, .55, .60) +
  draw_plot_label(
    c("","B"),
    c(.05, 0.55),
    c(1, 0.85),
    size = 16
  )

Removed 14 rows containing missing values (geom_point).

pgTop

pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c(“C”, ""), label_size = 14)

I don’t know whats going on below here

pgBottom <- pltDelta3  + geom_rect(aes(xmin = -2, xmax = -1.15, ymin = 170, ymax = 1000), colour = "gray50", fill = "white", inherit.aes = FALSE) + draw_plot(pltFlxLegend , -1.9, -575, .7)
pgBoth <- plot_grid(pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
                    pgBottom + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
                    ncol = 1, rel_heights = c(4, 4), labels = c("A", "C"), label_size = 16)

Removed 33 rows containing missing values (geom_point).

pgBoth

ggsave("figures/FluxDeepDive.png")

Saving 5 x 9 in image

ggsave("figures/FluxDeepDive.svg")

Saving 5 x 9 in image

Summary stats

Test for day to day and hourly variability in rate of change of flux (fifth root transformed)

There is variability with respect to depth, and day and hour Depth p = 0.03 R^2 = 0.088. Add affect of day p = 0.004, extra R^2 = 0.11, Add affect of hour p = 0.02 extra R2 = 0.12

bdsAddTime <- bds %>% 
  mutate(Hour = hour(time) + minute(time)/60, Day = day(time) + hour(time)/24 + minute(time)/24/60)

DFG1 <- gam(pracma::nthroot(DF/DZ, 5)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))

DFG2 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))

DFG3 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))

DFG_DayOnly <- gam(pracma::nthroot(DF/DZ, 5) ~  s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))

DFG_NoFifth <- gam(pracma::nthroot(DF/DZ, 1)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))

summary(DFG1)

Family: gaussian 
Link function: identity 

Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3) + 
    s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.18385    0.05998  -3.065  0.00407 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F p-value  
s(depth) 1.672  1.893 3.955  0.0606 .
s(Day)   1.907  1.991 4.661  0.0190 *
s(Hour)  0.774  2.000 0.658  0.1961  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.254   Deviance explained = 33.4%
GCV = 0.1732  Scale est. = 0.15112   n = 42

summary(DFG_DayOnly)

Family: gaussian 
Link function: identity 

Formula:
pracma::nthroot(DF/DZ, 5) ~ s(Day, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.18385    0.06485  -2.835  0.00722 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
         edf Ref.df     F p-value  
s(Day) 1.873  1.984 3.391  0.0404 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.129   Deviance explained = 16.8%
GCV = 0.18962  Scale est. = 0.17665   n = 42

summary(DFG_NoFifth)

Family: gaussian 
Link function: identity 

Formula:
pracma::nthroot(DF/DZ, 1) ~ s(depth, k = 3) + s(Day, k = 3) + 
    s(Hour, k = 4, bs = "cc")

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.028354   0.007924  -3.578 0.000979 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
               edf Ref.df     F p-value   
s(depth) 1.809e+00  1.963 9.247 0.00138 **
s(Day)   1.827e+00  1.970 3.361 0.06222 . 
s(Hour)  1.396e-08  2.000 0.000 0.71927   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.341   Deviance explained = 39.9%
GCV = 0.0029647  Scale est. = 0.0026375  n = 42

# summary(DFG2)
# summary(DFG3)
# 
# summary(DFG1)$r.sq - summary(DFG2)$r.sq
# summary(DFG2)$r.sq - summary(DFG3)$r.sq
# summary(DFG3)$r.sq

Time is now continuous day and hour

Plot of the gams above

#plot.new()
FluxGamPlot <- function(){
  par(mfrow = c(2,2))
  plot(DFG1)
  abline(h = 0, col = "gray30", lwd = 2)
  mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(1,1))
  abline(h = 0, col = "gray30", lwd = 2)
  mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(1,2))
  abline(h = 0, col = "gray30", lwd = 2)
  mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
}

FluxGamPlot()

png(filename = "./figures/FluxGamPlot.png", width = 10, height = 8, units = "in", res = 200)

FluxGamPlot()
dev.off()

png 
  2 

Figure 7

Difference from model expectations

(u mol C / m^3 / day)

disagFig <- bds %>% filter(project == "ETNP") %>%
  ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2) +
  scale_y_reverse(limits = c(1000, 0), breaks = ybreaks, labels = ylabels) +
  scale_x_continuous(limits = c(-1, 1)) +
  geom_vline(xintercept = 0) +   scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
  #labs(x = bquote("Observed - Modeled Small Particle Flux"~(μmol/m^3/day)), y = "Depth (m)") +
  labs(x = paste("Deviation from Model", expression((μmol/m^3/day)))) +
  scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  geom_hline(yintercept = PhoticBase, color = "darkgreen") +
  geom_hline(yintercept = OMZBase, color = "darkblue") +
  geom_hline(yintercept = OMZTop, color = "darkblue") #+
  #geom_hline(yintercept = DVMBase, color = "darkgoldenrod", lty = "dashed")
disagFig

#ggsave("..figures/FluxSizeShift.svg"

 ggsave("figures/FluxSizeShift.png")

Saving 6 x 4 in image

 ggsave("figures/FluxSizeShift.svg")

Saving 6 x 4 in image

Summary statistics

# bdsAddTime <- bds %>%
#   mutate(Hour = hour(time), Day = day(time))

bdsAddTime <- bds %>%
  mutate(Hour = hour(time) + minute(time)/60, Day = day(time) + hour(time)/24 + minute(time)/24/60)

OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= PhoticBase & depth <=500 & project == "ETNP"))

summary(OZG1)

Family: gaussian 
Link function: identity 

Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.08042    0.01049   7.665 1.33e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.584  1.827 19.185 4.77e-06 ***
s(Day)   1.771  1.946  8.926  0.00173 ** 
s(Hour)  1.266  2.000  1.964  0.05140 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.416   Deviance explained = 45.5%
GCV = 0.0082683  Scale est. = 0.0075948  n = 69

summary(OZG2)

Family: gaussian 
Link function: identity 

Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.08042    0.01079   7.451 2.85e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df      F  p-value    
s(depth) 1.558  1.805 18.393 8.14e-06 ***
s(Day)   1.838  1.974  6.924  0.00415 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.382   Deviance explained = 41.3%
GCV = 0.0085853  Scale est. = 0.0080382  n = 69

summary(OZG3)

Family: gaussian 
Link function: identity 

Formula:
ospsDZ ~ s(depth, k = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.08042    0.01169   6.877 2.59e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df     F  p-value    
s(depth) 1.453    1.7 17.31 2.23e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.274   Deviance explained =   29%
GCV = 0.0097831  Scale est. = 0.0094354  n = 69

summary(OZG1)$r.sq - summary(OZG2)$r.sq # Extra from Hour

[1] 0.03410746

summary(OZG2)$r.sq - summary(OZG3)$r.sq # Extrafrom Day

[1] 0.1074734

summary(OZG3)$r.sq # Depth

[1] 0.2741944

Plot of those gams Figure S10

OSMSGamPlot <- function(){
  par(mfrow = c(2,2))
  plot(OZG1)
  
  par(mfg = c(1,1))
  abline(h = 0, col = "gray30", lwd = 2)
  mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(1,2))
  abline(h = 0, col = "gray30", lwd = 2)
  mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
  par(mfg = c(2,1))
  abline(h = 0, col = "gray30", lwd = 2)
  mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
}

OSMSGamPlot()

png(filename = "./figures/OSMSGamPlot.png", width = 10, height = 8, units = "in", res = 200)

OSMSGamPlot()

dev.off()

png 
  2 

Figure 3

Trap data

trapFlux3 <- read_csv("dataOut/fluxMS_distilled.csv")

Parsed with column specification:
cols(
  Class = col_character(),
  Depth = col_double(),
  TrapID = col_character(),
  TrapType = col_character(),
  SampleType = col_character(),
  C_flux = col_double(),
  C_flux_umol = col_double()
)

UVPFluxComb <- read_csv("dataOut/CombinedProfileFluxEst_DS.csv")

Parsed with column specification:
cols(
  depth = col_double(),
  Flux = col_double()
)

UVPFluxOE <- read_csv("dataOut/ObservedVsExpectedFlux.csv")

Parsed with column specification:
cols(
  depth = col_double(),
  tn_flux = col_double(),
  profile = col_character(),
  project = col_character(),
  time = col_character(),
  tot_flux2 = col_double()
)

fluxMS_distilled_toPlot <- trapFlux3 %>%
  mutate(SampleType = recode(SampleType, `plus.p` = "plus-particles", top = "top-collector"))

Remove traps where mass spec didn’t work correctly 2-17 150 1-12 73m 1-12 148 2-14 100 |(TrapID == “2-17” & Depth == 150)

fluxMS_distilled_toPlot2 <- fluxMS_distilled_toPlot %>%
 filter(!((TrapID == "1-12") | (TrapID == "2-14" & Depth == 100)|(TrapID == "2-17" & Depth == 150)))
#fluxMS_distilled_toPlot2

Traps where mass spec didn’t work.

UVPFluxPlot00 <- UVPFluxComb %>% 
  ggplot(aes(y = depth))  + scale_y_reverse(limits = c(1000, 0)) +
  scale_x_continuous(limits = c(0, 200)) +
  geom_point(aes(y = Depth, x = C_flux_umol, shape = TrapType, ID = TrapID),
             colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot2) +
  geom_line(aes(x = Flux), size = 1, color = "black") +
  geom_point(aes(x = -1, y = -1, size = "UVP Estimate")) + # dummy point for the legend
  geom_point(aes(x = tot_flux2), size = 3, shape = 21, color = "white", fill = "black", data = UVPFluxOE) +
scale_shape_manual(values = c(25, 22))+
  scale_size_manual(values = 1, name = "") +
  ylab("Depth (m)") +
  #xlab(expression(Flux µmolC/m^2/day)) +
  xlab(expression(paste("x axis ", ring(A)^2))) +
  xlab(expression(paste("Flux (µ mol C/", m^2, "/day)"))) +
  
  guides(fill = guide_legend(override.aes = list(shape = 21))) +
  theme_cowplot() + 
  theme(
        legend.position = c(0.5, 0.4),
        legend.box.background = element_rect(color = "black", size = 0.5),
        legend.margin = margin(-10, 5, 10, 5)
  ) 

Ignoring unknown aesthetics: ID

# UVPFluxPlot <- UVPFluxPlot00 +
#   geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)

UVPFluxPlot00

ggsave("figures/FittedFlux.png")

Saving 7.29 x 4.5 in image

ggsave("figures/FittedFlux.svg")

Saving 7.29 x 4.5 in image

Figure S2

Example particle size distribution

TPPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = TotalParticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "TotalParticles Observed (#)")

nnpPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = n_nparticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() +
  #geom_vline(xintercept = 1) + geom_vline(xintercept = 5) +
  labs(y = "Depth (m)", x = "Binsize and Volume Normalized Particles (#/L/mm)")

FitPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = nnp_smooth, xmin = nnp_lower, xmax = nnp_upper, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() +
  #geom_vline(xintercept = 1) + geom_vline(xintercept = 5) +
  labs(y = "Depth (m)", x = "Smoothed - Normalized Particles (#/L/mm)") + geom_errorbar(width = 10, alpha = 0.5)

npLegend <- get_legend(FitPlot + theme(legend.box.margin = margin(0, 0, 40, 200)) + labs(col = expression(log[e](Size (mm)))))

Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).

plot_grid(
  TPPlot + theme(legend.position = "none"),
  nnpPlot + theme(legend.position = "none"),
  npLegend ,
  FitPlot + theme(legend.position = "none"), 
  labels = c("A", "B", "", "C")
)

Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).

ggsave("figures/AllParticleSizes.svg")

Saving 10 x 6.18 in image

ggsave("figures/AllParticleSizes.png")

Saving 10 x 6.18 in image

Figure 6

Weber Bianchi Figs

Just station 043

We are smoothing the station 043 data, with respect to depth and time. We are using this station because it is the only one to extend past 2000m. We find that the other profiles seem to give strange values near the ends of the ranges of the gams.

SameGam <- gam(TotalParticles ~s(log(lb), log(depth)), offset = log(vol * binsize), family = nb(),
    data = bes %>% filter(project == "ETNP", depth <= 2000, profile == "stn_043")) # Looks good!

gam.check(SameGam)

Method: REML   Optimizer: outer newton
full convergence after 4 iterations.
Gradient range [1.180573e-06,5.441536e-06]
(score 2243.109 & scale 1).
Hessian positive definite, eigenvalue range [5.359006,135.6228].
Model rank =  30 / 30 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

                      k' edf k-index p-value    
s(log(lb),log(depth)) 29  21    0.61  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

besE <- bes %>% filter(project == "ETNP")

lb_new <- exp(seq(from = log(0.1), to = log(2.1), by = 0.05))
ub_new <- lead(lb_new)
binsize_new <- ub_new - lb_new

lbbs <- tibble(lb = lb_new, ub = ub_new, binsize = binsize_new)

Expanded <- expand_grid(lb = exp(seq(from = log(0.1), to = log(2), by = 0.05)), depth = seq(from = 20, to = 2000, by = 20), time = as.factor(unique(besE$time))) %>% left_join(lbbs, by = "lb")

Pred <- exp(predict(SameGam, Expanded))

ToPlot <- bind_cols(Expanded, nnparticles = Pred) %>% mutate(time = as.character(time)) %>% mutate(nparticles = nnparticles * binsize)

WBColorMap <- ToPlot %>% filter(lb <= 2) %>%
  ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c(name = expression(log[10](Particles/m^3/mm))) + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs(y = "Depth (m)", x = "Size (mm)")
WBColorMap

mbGam <- ToPlot %>% group_by(depth)  %>% nest() %>%
  mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>% 
  mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2])) 

Particle size distribution, smoothed over all stations

pWBPSD <- mbGam %>%
  ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse()  + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept =  OMZBase, color = "darkblue") + labs(y = "Depth (m)", x = "Particle Size Distribution Slope")
pWBPSD

Large Particle Biomass

PubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb < 0.5) %>% group_by(depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup()
photicBiomass <- PubDf %>% filter(depth <= 165, depth >= 155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
PubDf <- PubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBS <- PubDf %>% 
  ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1.2)) + geom_hline(yintercept = 160, color = "darkgreen") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = OMZBase, color = "darkblue") + labs( x = "Small particle mass (norm.)")  + theme(axis.title.y = element_blank())
pWBS

#Small Particle Biomass

LubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb >= 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth)  %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- LubDf %>% filter(depth <= 165, depth >=155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
LubDf <- LubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBL <- LubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 160, color = "darkgreen") + labs( x = "Large particle mass (norm.)") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = OMZBase, color = "darkblue") + theme(axis.title.y = element_blank())
pWBL

Combine the three lower pannels

WBFig5 <- plot_grid(pWBPSD, pWBS,pWBL, nrow = 1, labels = c("B", "C", "D"))

Removed 6 row(s) containing missing values (geom_path).Removed 7 row(s) containing missing values (geom_path).

WBFig5

Four panel figure of Weber and Bianchi equivalent data

WBcombined <- plot_grid(WBColorMap + theme(plot.margin = unit(c(0,3,0, 3), "cm")), WBFig5, ncol = 1, labels = c("A", ""))
WBcombined

ggsave("figures/WBModelValidation.png")

Saving 8.5 x 6 in image

Figure S8

P16 Flux

Flux

scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxP16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = Flux_Smooth, group = factor(time)))  + geom_point(size = 3, stroke = 1)+
  geom_path() +
  scale_y_reverse(limits = c(1000, 0))+
  scale_x_log10(limits = c(35, 150),breaks = seq(from = 20, to = 150, by = 20)) +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + 
  scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
  geom_hline(yintercept = 200, color = "darkgreen") +
  theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm"))
# 
# 
# 
# pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
# pltFlxLegend <- get_legend(pltFlx)
# 
pltFlxP16

# #plotly::ggplotly(plt1)

Rate of change of flux – fifth root transformed

cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3P16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
  ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), group = factor(time)))  + geom_point(size = 3, stroke = 1)+
  geom_path() +
  scale_y_reverse(limits = c(1000, 0))+
  scale_x_continuous(limits = c(-1, .1), breaks = seq(from = -2, to = .75, by = 0.5)) +
  #scale_x_log10() +
   scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + 
  scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 200, color = "darkgreen")+
  labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")
  #labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")

pltDelta3P16

#plotly::ggplotly(plt1pos)

Difference from model

osms_p16 <- bds %>% filter(project == "P16") %>%
  ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), group = factor(time))) + geom_point(size = 3) + geom_path() + scale_y_reverse(limits = c(1000, 0)) +
  scale_x_continuous(limits = c(-1, 1)) +
  geom_vline(xintercept = 0) +   scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = "Observed - Modeled Small Particle Flux \n µmol/m^3/day") +
  scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = PhoticBase, color = "darkgreen") 
plotly::ggplotly(osms_p16)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

#ggsave("..figures/FluxSizeShift.svg"

Combine everything together

plot_grid(
  pltFlxP16,
  pltDelta3P16,
  osms_p16
)

Removed 28 rows containing missing values (geom_point).Removed 28 row(s) containing missing values (geom_path).Removed 32 rows containing missing values (geom_point).Removed 32 row(s) containing missing values (geom_path).Removed 29 rows containing missing values (geom_point).Removed 29 row(s) containing missing values (geom_path).

ggsave("figures/P16FluxRelate.svg")

Saving 8 x 8 in image

ggsave("figures/P16FluxRelate.png")

Saving 8 x 8 in image

Figure S10

Flux attenuation example

Take one profile, attenuate it, and show what it looks like

source("ModelStuff.R")

Parsed with column specification:
cols(
  Cruise = col_character(),
  Station = col_character(),
  `mon/dd/yyyy` = col_character(),
  `hh:mm` = col_time(format = ""),
  `Longitude [degrees east]` = col_double(),
  `Latitude [degrees north]` = col_double(),
  `Bottom Depth [m]` = col_double(),
  `Pressure [db]` = col_double(),
  `Temperature [degrees C]` = col_double(),
  `Temperature 2 [degrees C]` = col_double(),
  `Salinity [psu]` = col_double(),
  `Salinity 2 [psu]` = col_double(),
  `Fluorescense [mg/m^3]` = col_double(),
  `Beam Transmission [%]` = col_double(),
  PAR = col_double(),
  `Oxygen [umol/kg]` = col_double(),
  `Oxygen [% saturation]` = col_double()
)

scan_for_example <- bds %>% filter(project == "ETNP", depth < 500, depth > 200) %>% select(profile, depth, DFP, use_this_DFP, ospsDZ)

#loc_station = "stn_036"
loc_station = "stn_043"
loc_depth = 225
loc_prev_depth = 112.5

allDFPs <- bds %>% filter(profile == loc_station, depth >= loc_prev_depth, depth <= loc_depth) %>% summarize(DFP = prod(DFP), use_this_DFP = prod(use_this_DFP))

loc_DFP <-  allDFPs %>% pull(DFP)
loc_use_DFP <- allDFPs %>% pull(use_this_DFP)


for_single_disag <- bes %>% filter(profile == loc_station, depth %in% c(loc_prev_depth, loc_depth)) %>% select(depth, lb, nnp_smooth) %>%
  mutate(depth = recode(depth, `112.5` = "Shallow", `225` = "Deep")) %>% # I have no idea how to not hard code this bit
  pivot_wider(names_from = depth, values_from = nnp_smooth) 

with_disag <- for_single_disag %>%
  mutate(Predicted_Deep = remin_smooth_shuffle(Shallow, loc_use_DFP)) 
#remin_smooth_shuffle(for_single_disag$Shallow,loc_use_DFP)

for_plot_disag <- with_disag %>% pivot_longer(cols = -lb) %>% #filter(lb <= 5) %>%
  mutate(name = factor(name, levels = c("Shallow", "Deep", "Predicted_Deep"))) %>%
  mutate(name = recode_factor(name, Shallow = "Shallow (112.5m)", Deep = "Deep (225m)", Predicted_Deep = "Predicted Deep (225m)"))

for_plot_disag %>% ggplot(aes(x = lb, y = value, shape = name)) + geom_point() + scale_x_log10() + scale_y_log10() + scale_shape_manual(values = c(1, 6, 3)) + theme(legend.title = element_blank()) + labs(x = "Particle Size (mm)", y = "Normalized Particle Abundance (#/L/mm)")

ggsave("figures/DisagExample.png")

Saving 7.29 x 4.5 in image

ggsave("figures/DisagExample.svg")

Saving 7.29 x 4.5 in image